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ABSTRACT 

We describe results from a fully self-consistent three dimensional hydrodynamical simulation of the 
formation of one of the first stars in the Universe. Dark matter dominated pre-galactic objects form 
because of gravitational instability from small initidal density perturbations. As they assemble via 
hierarchical merging, primordial gas cools through ro-vibrational lines of hydrogen molecules and sinks 
to the center of the dark matter potential well. The high redshift analog of a molecular cloud is formed. 
When the dense, central parts of the cold gas cloud become self-gravitating, a dense core of ~ 100 Mq 
undergoes rapid contraction. At densities n > 10° cm ~ 3 a 1 Mq proto-stellar core becomes fully molecular 
due to three-body H 2 formation. Contrary to analytical expectations this process does not lead to 
renewed fragmentation and only one star is formed. The calculation is stopped when optical depth 
effects become important, leaving the final mass of the fully formed star somewhat uncertain. At this 
stage the protostar is acreting material very rapidly (~ lO _ 2 M 0 yr _1 ). Radiative feedback from the 
star will not only halt its growth but also inhibit the formation of other stars in the same pre-galactic 
object (at least until the first star ends its life, presumably as a supernova). We conclude that at most 
one massive (M 1M@) metal free star forms per pre-galactic halo, consistent with recent abundance 
measurements of metal poor galactic halo stars. 

1. MOTIVATION 

Chemical elements heavier than Lithium are synthesized 
in stars. Such “metals” are observed at times when the 
Universe was only < 10% of its current age in the inter 
galactic medium (IGM) as absorption lines in quasar spec¬ 
tra (see Ellison et al. 2000, and references therein). Hence, 
these heavy elements not only had to be synthesized but 
also released and distributed in the IGM within the first 
billion years. Only supernovae of sufficiently short lived 
massive stars are known to provide such an enrichment 
mechanism. This leads to the prediction that 

the first generation of cosmic structures formed massive 
stars (although not necessarily only massive stars). 

In the past 30 years it has been argued that the first 
cosmological objects form globular clusters ( 1 ), super 
massive black holes (2), or even low mass stars ( 3 ). This 
disagreement of theoretical studies might at first seem sur¬ 
prising. However, the first objects form via the gravita¬ 
tional collapse of a thermally unstable reactive medium, 
inhibiting conclusive analytical calculations. The prob¬ 
lem is particularly acute because the evolution of all other 
cosmological objects (and in particular the larger galaxies 
that follow) will depend on the evolution of the first stars. 

Nevertheless, in comparison to present day star forma¬ 
tion, the physics of the formation of the first star in the 
universe is rather simple. In particular: 

• the chemical and radiative of processes in the 
primordial gas are readily understood. 


• strong magnetic fields are not expected to exist at 
early times. 

• by definition no other stars exist to influence the 
environment through radiation, winds, supernovae, 
etc. 

• the emerging standard model for structure 
formation provides appropriate initial conditions. 

In previous work we have presented three-dimensional 
cosmological simulations of the formation of the first ob¬ 
jects in the universe (4, 5) including first applications of 
adaptive mesh refinement (AMR) cosmological hydrody¬ 
namical simulations to first structure formation (6, 7, ABN 
hereafter) . In these studies we achieved a dynamic range 
of up to 2 x 10 5 and could follow in detail the formation of 
the first dense cooling region far within a pre-galactic ob¬ 
ject that formed self-consistently from linear density fluc¬ 
tuation in a cold dark matter cosmology. Here we report 
results from simulations that extend our previous work by 
another 5 orders of magnitude in dynamic range. For the 
first time it is possible to bridge the wide range between 
cosmological and stellar scale. 

2. SIMULATION SETUP AND NUMERICAL ISSUES 

We employ an Eulerian structured adaptive mesh re¬ 
finement cosmological hydrodynamical code developed by 
Bryan and Norman (9, 10). The hydrodynamical equa¬ 
tions are solved with the second order accurate piecewise 
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parabolic method ( 11 ; 12) where a Riemann solver en¬ 
sures accurate shock capturing with a minimum of numer¬ 
ical viscosity. We use initial conditions appropriate for a 
spatially flat Cold Dark Matter cosmology with 6% of the 
matter density contributed by baryons, zero cosmological 
constant, and a Hubble constant of 50 km/s/Mpc (8). The 
power spectrum of initial density fluctuations in the dark 
matter and the gas are taken from the computation by the 
publicly available Boltzmann code CMBFAST (13) at red- 
shift 100 (assuming an Harrison-Zel’dovich scale—invariant 
initial spectrum). 

We set up a three dimensional volume with 128 comov¬ 
ing kpc on a side and solve the cosmological hydrodynam¬ 
ics equations assuming periodic boundary conditions. This 
small volume is adequate for our purpose, because we are 
interested in the evolution of the first pre-galactic object 
within which a star may be formed by a redshift of z ~ 20. 
First we identify the Lagrangian volume of the first proto- 
galactic halo with a mass of ~ 10 6 M 0 in a low resolution 
pure N-body simulation. Then we generate new initial 
conditions with four initial static grids that cover this Lan- 
grangian region with progressively finer resolution. With a 
64 3 top grid and a refinement factor of 2 this specifies the 
initial conditions in the region of interest equivalent to a 
512 3 uni-grid calculation. For the adopted cosmology this 
gives a mass resolution of 1.1 M 0 for the dark matter (DM, 
hereafter) and 0.07 M 0 for the gas. The small DM masses 
ensure that the cosmological Jeans mass is resolved by at 
least ten thousand particles at all times. Smaller scale 
structures in the dark matter will not be able to influence 
the baryons because of their shallow potential wells. The 
theoretical expectation holds, because the simulations of 
ABN which had 8 times poorer DM resolution led to iden¬ 
tical results on large scales as the simulation presented 
here. 

During the evolution, refined grids are introduced with 
twice the spatial resolution of the parent (coarser) grid. 
These child (finer) meshes are added whenever one of three 
refinement criteria are met. Two Langrangian criteria en¬ 
sure that the grid is refined whenever the gas (DM) density 
exceeds 4.6 (9.2) its initial density. Additionally, the local 
Jeans length is always covered by at least 64 grid cells 1 
(4 cells per Jeans length would be sufficient, 14)- We have 
also carried out the simulations with identical initial con¬ 
ditions but varying the refinement criteria. In one series 
of runs we varied the number of mesh points per Jeans 
length. Runs with 4, 16, and 64 zones per Jeans length 
are indistinguishable in all mass weighted radial profiles 
of physical quantities. No change in the angular momen¬ 
tum profiles could be found, suggesting negligible numer¬ 
ical viscosity effects on angular momentum transport. A 
further refinement criterion that ensured the local cooling 
time scale to be longer than the local Courant time also 
gave identical results. This latter test checked that any 
thermally unstable region was identified. 

The simulation follows the non-equilibrium chemistry 
of the dominant nine species species (H, H + , H - , e - , 
He, He+, He++, H 2 , and H+) in primordial gas. Fur¬ 
thermore, the radiative losses from atomic and molecular 
line cooling, Compton cooling and heating of free elec¬ 
trons by the cosmic background radiation are appropri¬ 


ately treated in the optically thin limit (15, 16). To ex¬ 
tend our previous the studies to higher densities three es¬ 
sential modifications to the code were made . First we 
implemented the three-body molecular hydrogen forma¬ 
tion process in the chemical rate equations. For tem¬ 
peratures below 300 K we fit to the data of Orel (17) 
to get k 3 b = 1-3 x 10 _32 (T/300 K) -0 - 38 cm 6 s _1 . Above 
300 K we then match it continuously to a powerlaw (3) 
k^b = 1.3xl0 -32 (T/300 K) -1 cm 6 s -1 . Secondly, we intro¬ 
duce a variable adiabatic index for the gas (18). The dissi¬ 
pative component (baryons) may collapse to much higher 
densities than the collisionless component (DM). The dis¬ 
crete sampling of the DM potential by particles can then 
become inadequate and result in artificial heating of the 
baryons (cooling for the DM) once the gas density becomes 
much larger than the local DM density. To avoid this, we 
smooth the DM particles with a Gaussian of width 0.05 pc 
for grids with cells smaller than this length. At this scale, 
the enclosed gas mass substantially exceeds the enclosed 
DM mass. 

The standard message passing library (MPI) was used to 
implement domain decomposition on the individual levels 
of the grid hierarchy as a parallelization strategy. The code 
was run in parallel on 16 processors of the SGI Origin2000 
supercomputer at the National Center for Supercomput¬ 
ing Applications at the University of Illinois at Urbana 
Champaign. 

We stop the simulation at a time when the molecular 
cooling lines reach an optical depth of ten at line cen¬ 
ter because our numerical method cannot treat the diffi¬ 
cult problem of time-dependent radiative line transfer in 
multi-dimensions. At this time the code utilizes above 
5500 grids on 27 refinement levels with 1.8 x 10 7 « 260 3 
computational grid cells. An average grid therefore con¬ 
tains ~ 15 3 cells. 

3. RESULTS 

3.1. Characteristic mass scales 

Our simulations (Fig. 1, Fig. 2), identify at least four 
characterise mass scales. From the outside going in, one 
observes infall and accretion onto the pre-galactic halo 
with a total mass of 7 x 10 5 M 0 , consistent with previous 
studies (4, 19, 5, ABN, and 20 for discussion and refer¬ 
ences). 

At a mass scale of about 4000 solar mass (r ~ 10 pc) 
rapid cooling and infall is observed. This is accompanyed 
by the first of three valleys in the radial velocity distribu¬ 
tion (Fig. 2E). The temperature drops and the molecular 
hydrogen fraction increases. It is here, at number densities 
of ~ 10 cm -3 , that the high redshift analog of a molecular 
cloud is formed. Although the molecular mass fraction is 
not even 0.1% it is sufficient to cool the gas rapidly down 
to ~ 200 K. The gas cannot cool below this temperature 
because of the sharp decrease in the cooling rate below 
~ 200AT. 

At redshift 19 (Fig. 2), there are only two mass scales; 
however, as time passes the central density grows and even¬ 
tually passes 10 4 cm -3 , at which point the ro-vibrational 
levels of H 2 are populated at their equilibrium values and 
the cooling time becomes independent of density (instead 
of inversely proportional to it). This reduced cooling effi- 


The Jeans mass which is the relevant mass scale for collapse and fragmentation is thus resolved by at least 47r32 3 /3 rs 1.4 X 10® cells. 
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Fig. 1.— Overview of the evolution and collapse forming a primordial star in the universe. The top row shows projections of the gas 
density of one thousands of the simulation volume approximately centered at the pre-galactic object within which the star is formed. The 
four projections from left to right are taken at redshifts 100, 24, 20.4, and 18.2 respectively. The pre—galactic objects form from very small 
density fluctuations, and continously merge to form larger objects. The middle and bottom row show thin slices through the gas density and 
temperature at the final simulation output. The leftmost panels are on the scale of the simulation volume ~ 6 proper kpc. The panels to the 
right zoom in towards the forming star and have side lengths of 600 pc, 6pc, and 0.06 pc (12000 astronomical units). The color maps (going 
from black to blue, green, red, yellow) are logarithmic and the associated values were adjusted considerably to visualize the ~ 17 orders of 
magnitude in density covered by these simulations. In the left panels the larger scale structures of filaments and sheets are seen. At their 
intersections a pre-galactic object of ~ 10 6 Mq is formed. In the temperature slice (second panel - bottom row) one sees how the gas shock 
heats as it falls into the pre-galactic object. After passing the acretion shock the material forms hydrogen molecules and starts to cool. The 
cooling material accumulates at the center of the object and forms the high redshift molecular cloud analog (third panel from the right) which 
is dense and cold (T ~ 200 K). Deep within the molecular cloud a few hundred Kelvin warmer core of ~ 100 Mq is formed (right panel) 
within which a 1 Mq proto-star is formed (yellow region in the right panel of the middle row). 


ciency leads to an increase in the temperature (Fig. 2D). 
As the temperature rises, the cooling rate again increases 
(it is 1000 times higher at 800 K than at 200 K), and the 
inflow velocities slowly climb. 

In order to better understand what happens next, we 
examine the stability of an isothermal gas sphere. The 
critical mass for gravitational collapse given an external 
pressure P ex t (BE mass hereafter) is given by Ebert {21) 
and Bonnor {22) as: 


Mbe = 1.18 Mq 


d p —1/2. 

Q 3/2 ext ) 


dP_ 

dp 


lk B T 

prriH 


(1) 


Here P ex t is the external pressure and G, k B , and c s the 
gravitational constant, the Boltzmann constant and the 
sound speed, respectively. We can estimate this critical 


mass locally if we set the external pressure to be the lo¬ 
cal pressure to find Mbe ~ 2 OM 0 T 3 / 2 n _ 1 / 2 /i _ 2 7 2 where 
p ss 1.22 is the mean mass per particle in units of the pro¬ 
ton mass. Using an adiabatic index 7 = 5/3, we plot the 
ratio of the enclosed gas mass to this modified BE mass in 
Figure 3. 

Our modeling shows (Fig. 3), that by the fourth con¬ 
sidered output time, the central 100 M 0 exceeds the BE 
mass at that radius, indicating unstable collapse. This is 
the third mass scale and corresponds to the second local 
minimum in the radial velocity curves (Fig. 2E). The in¬ 
flow velocity is 1 km s ” 1 is still subsonic. Although this 
mass scale is unstable, it does not represent the smallest 
scale of collapse in our simulation. This is due to the in- 
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Fig. 2.— Radial mass weighted averages of various physical quantities at seven different output times. Panel A shows the evolution of the 
particle number density in cm -3 as a function of radius at redshift 19 (solid line), nine Myrs later (dotted lines with circles), 0.3Myr later 
(dashed line), 3 x 10 4 yr later (long dashed line), 3 x 10 3 yr later (dot-dashed line), 1.5 x 10 3 yr later (solid line) and finally 200 years later 
(dotted line with circles) at z = 18.181164. The two lines between 10 -2 and 200 pc give the DM mass density in GeV cm -3 at z=19 and the 
final time, respectively. Panel B gives the enclosed gas mass as a function of radius. In C the mass fractions of atomic hydrogen and molecular 
hydrogen are shown. Panel D and E illustrate the temperature evolution and the mass weighted radial velocity of the baryons, respectively. 
The bottom line with filled symbols in panel E shows the negative value of the local speed of sound at the final time. In all panels the same 
output times correspond to the same line styles. 


creasing molecular hydrogen fraction. 

When the gas density becomes sufficiently large (~ 
10 10 cm -3 ), three-body molecular hydrogen formation be¬ 
comes important. This rapidly increases the molecu¬ 
lar fraction (Fig. 2C) and hence the cooling rate. The 
increased cooling leads to lower temperatures and even 
stronger inflow and. At a mass scale of ~ 1M 0 , not only 
is the gas nearly completely molecular, but the radial in¬ 
flow has become supersonic (Fig. 2E). When the H 2 mass 
fraction approaches unity, the increase in the cooling rate 
saturates, and the gas goes through a radiative shock. This 
marks the first appearance of the proto-stellar accretion 
shock at a radius of about 20 astronomical units from its 
center. 


3.2. Chemo-Thermal Instability 

When the cooling time becomes independent of density 
the classical criterion for fragmentation t coo i < td yn oc 
n -1 / 2 (23) cannot be satisfied at high densities. However, 
in principal the medium may still be subject to thermal 


instability. The instability criterion is 


P 


dL\ 

3P J T=const. 



p=const. 


+ L(p, T) > 0, (2) 


where L denotes the cooling losses per second of a fluid 
parcel and T and p are the gas temperature and mass 
density, respectively. At densities above the critical densi¬ 
ties of molecular hydrogen the cooling time is independent 
of density, i.e. dL/dp = A(T) where A (T) is the high den¬ 
sity cooling function (e.g. 24). Fitting the cooling function 
with a power-law locally around a temperature To so that 
A (T) oc (T/T 0 ) a one finds dL/dT = paA(T)/T. Hence, 
under these circumstances the medium is thermally stable 
if a > 2. Because, a > 4 for the densities and tempera¬ 
tures of interest, we conclude that the medium is thermally 
stable. The above analysis neglects the heating from con¬ 
traction, but this only strengthens the conclusion. If heat¬ 
ing balances cooling one can neglect the +L(p , T) term in 
equation (2) and find the medium to be thermally stable 
for a > 1. 
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Fig. 3. — Ratio of enclosed gas mass to the locally estimated 
Bonnor—Ebert mass ( Mg e rs 61 M 0 T^ 2 n _ 1 / 2 /i _2 ) for various out¬ 
put times. The enclosed gas mass exceeds the BE mass at two dif¬ 
ferent mass scales, ~ 1M@ and ~ 100Mg. The line-styles in the 
Figure correspond to the output times shown in Fig. 2. 

However, here we neglected the chemical processes. The 
detailed analysis for the case when chemical processes oc¬ 
cur on the collapse time-scale is well known (25). This 
can be applied to primordial star formation (26) includ¬ 
ing the three-body formation of molecular hydrogen (3) 
which drives a chemo-thermal instability. Evaluating all 
the terms in this modified instability criterion (26, equa¬ 
tion 36) one finds the simple result that for molecular mass 
fractions / < 6/(2 a + 1) the medium is expected to be 
chemo-thermally unstable. These large molecular frac¬ 
tions illustrate that the strong density dependence of the 
three body H 2 formation dominates the instability. Ex¬ 
amining the three dimensional temperature and H 2 den¬ 
sity field we clearly see this chemo-thermal instability at 
work. Cooler regions have larger H 2 fractions. How¬ 
ever, no corresponding large density inhomogeneities are 
found and fragmentation does not occur. This happens 
because of the short sound crossing times in the collapsing 
core. When the H 2 formation time-scale becomes shorter 
than the cooling time the instability originates. However, 
as long as the sound crossing time is much shorter than 
the chemical and cooling time scales the cooler parts are 
efficiently mixed with the warmer material. This holds 
in our simulation until the final output where for the 
first time the H 2 formation time scale becomes shorter 
than the sound-crossing time. However, at this point the 
proto-stellar core is fully molecular and stable against the 
chemo-thermal instability. Consequently no large density 
contrasts are formed. Because at these high densities the 
optical depth of the cooling radiation becomes larger than 
unity the instability will be suppressed even further. 

3.3. Angular momentum 

Interestingly, rotational support does not halt the col¬ 
lapse. This is for two reaons. The first is shown in panel 
A of Fig. 4, which plots the specific angular momentum 
against enclosed mass for the same seven output times dis¬ 
cussed earlier. Concentrating on the first output (Fig. 4), 
we see that the central gas begins the collapse with a spe¬ 


cific angular momentum only ~ 0.1% as large as the mean 
value. This type of angular momentum profile is typical 
of halos produced by gravitational collapse (e.g. 27), and 
means that the protostellar gas starts out without much 
angular momentum to lose. As a graphic example of this, 
consider the central one solar mass of the collapsing region. 
It has only an order of magnitude less angular momentum 
at densities n > 10 13 cm -3 than it had at n > 10 6 cm -3 
although it collapsed by over a factor 100 in radius. 

The remaining output times (Fig. 4) indicate that there 
is some angular momentum transport within the central 
100 M 0 (since L plotted as a function of enclosed mass 
should stay constant as long as there is no shell crossing). 
In panel C, we divide L by r to get a typical rotational ve¬ 
locity and in panels B and D compare this velocity to the 
Keplerian rotational velocity and the local sound speed, 
respectively. 

We find that the typical rotational speed is a factor 
two to three below that required for rotational support. 
Furthermore, we see that this azimuthal speed never sig¬ 
nificantly exceeds the sound speed, although for most the 
mass below 100 M 0 it is comparable in value. We interpret 
this as evidence that it is shock waves during the turbulent 
collapse that are responsible for much of the transported 
angular momentum. A collapsing turbulent medium is 
different from a disk in Keplarian rotation. At any radius 
there will be both low and high angular momentum mate¬ 
rial, and pressure forces or shock waves can redistribute 
the angular momentum between fluid elements. Lower 
angular momentum material will selectively sink inwards, 
displacing higher angular momentum gas. This hydrody¬ 
namic transport of angular momentum will be suppressed 
in situation where the collapse proceeds on the dynamical 
time rather on the longer cooling time as in the presented 
case. This difference in cooling time and the widely dif¬ 
ferent initial conditions may explain why this mechanism 
has not been observed in simulations of present day star 
formation (e.g. 28, and references therein). However, such 
situations may also arise in the late stages of the formation 
of present day stars and in scenarios for the formation of 
super-massive black holes. 

To ensure that the angular momentum transport is not 
due to numerical shear viscosity (29) we have carried out 
the resolution study discussed above. We have varied the 
effective spatial resolution by a factor 16 and found iden¬ 
tical results. Furthermore, we have run the adaptive mesh 
refinement code with two different implementations of the 
hydrodynamics solver. The resolution study and the re¬ 
sults presented here were carried out with a direct piece- 
wise parabolic method adopted for cosmology (11; 12). 
We ran another simulation with the lower order ZEUS hy¬ 
drodynamics (30) and still found no relevant differences. 
These tests are not strict proof that the encountered an¬ 
gular momentum transport is not caused by numerical ef¬ 
fects; however, they are reassuring. 

3.4. Magnetic Fields? 

The strength of magnetic fields generated around the 
epoch of recombination is minute. In contrast, phase 
transitions at the qantum-chromo-dynamic (QCD) and 
electro-weak scales may form even dynamically impor¬ 
tant fields. While there is a plethora of such scenarios 
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Fig. 4.— Radial mass weighted averages of various physical quantities related to the angular momentum of the gas. The seven different 
output times correspond to the ones described in Fig. 2. Panel A shows the specific angular momentum L in km s —1 pc as a function of 
enclosed gas mass. The typical rotational speed L/r is shown in panel C and its ratio to the Keplerian velocity Vpcepler — (GM/r) 1 / 2 and 
the local speed of sound in panel B and D, respectively. 


for primordial magnetic field generation in the early uni¬ 
verse they are not considered to be an integral part of our 
standard picture of structure formation. This is because 
not even the order of these phase transitions is known 
(31), and references therein). Unfortunately, strong pri¬ 
mordial small-scale (<C 1 comoving Mpc) magnetic fields 
are poorly constrained observationally (32). 

The critical magnetic field for support of a cloud (33) 
allows a rough estimate up to which primordial magnetic 
field strengths we may expect our simulation results to 
hold. For this we also assume a flux frozen flow with no 
additional amplification of the magnetic field other than 
the contraction (B oc p 2 / 3 ). For a comoving B field of 
> 3 x 10“ 11 G on scales < 100 kpc the critical field needed 
for support may be reached during the collapse possibly 
modifying the mass scales found in our purely hydrody¬ 
namic simulations. However, the ionized fraction drops 
rapidly during the collapse because of the absence of cos¬ 
mic rays ionizations. Consequently ambipolar diffusion 
should be much more effective in the formation of the first 
stars even if such strong primordial magnetic fields were 
present. 


4. DISCUSSION 

Previously we discussed the formation of the pre- 
galactic object and the primordial “molecular cloud” that 
hosts the formation of the first star in the simulated patch 
of the universe (7). These simulations had a dynamic 
range of ~ 10 5 and identified a ~ 100 M 0 core within the 
primordial “molecular cloud” undergoing renewed gravi¬ 
tational collapse. The fate of this core was unclear be¬ 
cause there was the potential caveat that three body H 2 
formation could have caused fragmentation. Indeed this 
further fragmentation had been suggest by analytic work 
(26) and single zone models (3). The three dimensional 
simulations described here were designed to be able to test 
whether the three body process will lead to a break up of 
the core. No fragmentation due to three body H 2 forma¬ 
tion is found. This is to a large part because of the slow 
quasi-hydrostatic contraction found in ABN which allows 
sub-sonic damping of density perturbations and yields a 
smooth distribution at the time when three body H 2 for¬ 
mation becomes important. Instead of fragmentation a 
single fully molecular proto-star of ~ 1M 0 is formed at 
the center of the ~ 100 M 0 core. 

However, even with extraordinary resolution, the final 
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mass of the first star remains unclear. Whether all the 
available cooled material of the surroundings will accrete 
onto the proto-star or feedback from the forming star will 
limit the further accretion and hence its own growth is 
difficult to compute in detail. 

Within 10 4 yr about 70 M 0 may be accreted assuming 
that angular momentum will not slow the collapse (Fig. 5). 
The maximum of the accretion time of ~ 5 x 10 6 yr is 
at ~ 600 Mg. However, stars larger than 100 M 0 will 
explode within ~ 2 Myr. Therefore, it seems unlikely 
(even in the absence of angular momentum) that there 
would be sufficient time to accrete such large masses. 



Fig. 5.— The accretion time as function of enclosed gas mass. 
The line with symbols gives M(r)/[4rrp(r)r 2 |u r (r)|]. The solid line 
simply shows how long it would take the mass to move to r = 0 if 
it were to to keep its current radial velocity ( r/v r {M )). 

A one solar mass proto-star will evolve too slowly to 
halt substantial accretion. From the accretion time pro¬ 
file (Fig. 5) one may argue that a more realistic minimum 
mass limit of the first star should be >, 30 M 0 because this 
amount would be accreted within a few thousand years. 
This is a very short time compared to expected proto- 
stellar evolution times. However, some properties of the 
primordial gas may make it easier to halt the accretion. 
One possibility is the destruction of the cooling agent, 
molecular hydrogen, without which the acreting material 
may reach hydrostatic equilibrium. This may or may not 
be sufficient to halt the accretion. One may also imag¬ 
ine that the central material heats up to 10 4 K, allowing 
Lyman-a cooling from neutral hydrogen. That cooling re¬ 
gion may expand rapidly as the internal pressure drops 
because of infall, possibly allowing the envelope to accrete 
even without molecular hydrogen as cooling agent. Addi¬ 
tionally, radiation pressure from ionizing photons as well 
as atomic hydrogen Lyman series photons may become 
significant and eventually reverse the flow. The mecha¬ 
nisms discussed by Haehnelt (1993) on galactic scales will 
play an important role for the continued accretion onto the 
proto-star. This is an interplay of many complex physical 
processes because one has a hot ionized Stromgren sphere 
through which cool and dense material is trying to accrete. 
In such a situation one expects a Raleigh-Taylor type in¬ 
stability that is modified via the geometry of the radiation 


held. 

At the final output time presented here there are ~ 
4 x 10 57 hydrogen molecules in the entire protogalaxy. 
Also the H 2 formation time scale is long because there 
are no dust grains and the free electrons (needed as a 
catalyst) have almost fully recombined. Hence, as soon 
as the the first UV photons of Lyman Werner band fre¬ 
quencies are produced there will be a rapidly expanding 
photo-dissociating region (PDR) inhibiting further cool¬ 
ing within it. This photo-dissociation will prevent fur¬ 
ther fragmentation at the molecular cloud scale. I.e. no 
other star can be formed within the same halo before the 
first star dies in a supernova. The latter, however, may 
have sufficient energy to unbind the entire gas content of 
the small pre-galactic object it formed in (35). This may 
have interesting feedback consequences for the dispersal of 
metals, entropy and magnetic Held into the intergalactic 
medium (36, 37). 

Smoothed particle hydrodynamics (SPH, e.g. 38), used 
extensively in cosmological hydrodynamics, has been em¬ 
ployed (20) to follow the collapse of solid body rotating 
uniform spheres. The assumption of coherent rotation 
causes these clouds to collapse into a disk which devel- 
opes filamentary structures which eventually fragment to 
form dense clumps of masses between 100 and 1000 solar 
masses. It has been argued that these clumps will con¬ 
tinue to accrete and merge and eventually form very mas¬ 
sive stars. These SPH simulation have unrealistic initial 
conditions and much less resolution then our calculations. 
However, they also show that many details of the collapse 
forming a primordial star are determined by the properties 
of the hydrogen molecule. 

We have also simulated different initial density fields 
for a Lambda CDM cosmology. There we have focused 
on halos with different clustering environments. Although 
we have not followed the collapse in these halos to proto- 
stellar densities, we have found no qualitative differences 
in the “primordial molecular cloud” formation process as 
discussed in ABN. Also other AMR simulations (39) give 
consistent results on scales larger than 1 pc. In all cases 
a cooling flow forms the primordial molecular cloud at 
the center of the dark matter halo. We conclude that 
the molecular cloud formation process seems to be inde¬ 
pendent of the halo clustering properties and the adopted 
CDM type cosmology. Also the mass scales for the core 
and the proto-star are determined by the local Bonnor 
Ebert mass. Consequently, we expect the key results dis¬ 
cussed here to be insensitive to variations in cosmology or 
halo clustering. 

5. CONCLUSION 

The picture arising from these numerical simulations has 
some very interesting implications. It is possible that all 
metal free stars are massive and form in isolation. Their 
supernovae may provide the metals seen in even the low¬ 
est column density quasar absorption lines ( 40 , and ref¬ 
erences therein). Massive primordial stars offer a natural 
explanation for the absence of purely metal free low mass 
stars in the Milky Way. The consequences for the forma¬ 
tion of galaxies may be even more profound in that the 
supernovae provide metals, entropy, and magnetic fields 
and may even alter the initial power spectrum of density 







fluctuations of the baryons. 

Interestingly, it has been recently argued, from abun¬ 
dance patterns, that in low metallicity galactic halo stars 
seem to have been enriched by only one population of mas¬ 
sive stars (41)- These results, if confirmed, would represent 
strong support for the picture arising from our ab initio 
simulations of first structure formation. 

To end on a speculative note there is suggestive evi¬ 
dence that links gamma ray bursts to sites of massive star 
formation (e.g. 4%)- It would be very fortunate if a sig¬ 
nificant fraction of the massive stars naturally formed in 
the simulations would cause gamma ray bursts (e.g. ^3). 
Such high redshift bursts would open a remarkably bright 
window for the study of the otherwise dark (faint) Ages. 
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